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Numerical solution of multi-order 
fractional differential equations via the 
sinc collocation method 


E. Hesameddini* and E. Asadollahifard 


Abstract 


In this paper, the sinc collocation method is proposed for solving linear 
and nonlinear multi-order fractional differential equations based on the new 
definition of fractional derivative which is recently presented by Khalil, R., 
Al Horani, M., Yousef, A. and Sababeh, M. in A new definition of fractional 
derivative, J. Comput. Appl. Math. 264 (2014), 65-70. The properties 
of sinc functions are used to reduce the fractional differential equation to a 
system of algebraic equations. Several numerical examples are provided to 
illustrate the accuracy and effectiveness of the presented method. 


Keywords: Sinc function; Fractional differential equations; Multi-order 
FDEs; Collocation method. 


1 Introduction 


One of the old fields of mathematics is fractional calculus which dates back 
to the time of Leibniz [1] and from then many studies were done in this field 
[14]-[12]. Fractional differential equations (FDEs) have attracted the interest 
of researchers in many areas such as Physics, Chemistry, Engineering and 
Social Sciences [22, 15]. The analytic results on the existence and uniqueness 
of solutions to the FDEs have been investigated by many authors [11, 22, 16]. 
Generally, most of the FDEs do not have analytic solutions, so one has to 
resort to approximation and numerical methods. 

One class of FDEs is multi-order fractional differential equations. They 
have been used to model various types of visco-elastic damping [22] and are 
expressed as follows 
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DO y(x) = F(x, y(x), D® y(a), ..., D& y(x)), «ce I=[0,l], (1) 

with initial conditions 
D}(0)=d;, i=0,1,..,m—1, meN (2) 


where m—1<a<m, 0 < $, < Bo <... < Be < q@ and the values of 
d; (i = 0,1,...,m—1) describe the initial state of y(x). D‘“y indicates the 
fractional derivative of order a of y. Up to now, whenever this equation 
was under study, in most cases the fractional derivative was in the sense of 
Caputo definition. In this paper, we imply the new definition of conformable 
fractional derivative [18] which will be defined later. Depending on F’,, this 
equation classifies into linear and nonlinear. 


In [14], it has been proved that equation (1) subject to the initial condi- 
tions (2) and under natural Lipschitz conditions imposed on F' has a unique 
continuous solution. 


Since the last decade, extensive research has been conducted on the devel- 
opment of numerical methods for equation (1). Doha et al.[25] proposed an 
efficient spectral tau and collocation method based on the Chebyshev poly- 
nomials for solving this equation. Extension of the tau method based on the 
shifted Legendre Gauss-Lobbato quadrature is used for solving equation (1) 
in [9]. In [12], this equation is converted into a system of FDEs and the 
shifted Chebyshev operational matrix method is used to solve the resultant 
system. Some other works on this problem are: piecewise polynomial col- 
location [17], Haar wavelet method [20], Lagrange wavelet method [23] and 
second kind Chebyshev wavelet method [30]. 


In this work, we apply the sinc collocation method for solving equation 
(1). The sinc method is an efficient method developed by Stenger [24]. It was 
widely used for the numerical solution of initial and boundary value problems 
[13, 19, 8], not only because of its exponential convergence rate but also due 
to its ability in handling problems with singularities. To the best of our 
knowledge, the sinc collocation method has not been used for solving FDEs 
directly. In this work, based on the new definition of fractional derivative 
[18], we compute the fractional derivative of the sinc function and apply it 
for solving equation (1). 


The remainder of this paper is organized as follows: in Section 2, some 
definitions and theorems are presented that will be used in later sections. The 
proposed method is discussed in Section 3. Section 4 is devoted to numerical 
experiments. Finally some remarks are concluded. 
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2 Preliminaries 


In this section, we recall some necessary definitions and mathematical pre- 
liminaries of the fractional theory and sinc method which will be used further 
in this paper. 


2.1. The fractional derivative 

The fractional calculus involves different definitions of fractional derivative 
operators such as Caputo and Riemmann-Lioville fractional derivative[22, 1]. 
One of the most recent works on the theory of derivatives of fractional order 
is done by Khalil et al. [18] which is the simplest definition. Up to now, some 
works were done based on this new definition [1, 2, 22]. In what follows, at 
first the conformable fractional derivative is defined and then some fantastic 
properties of this definition are presented. 


Definition 1. [18] Let a € (n,n +1], and f be an n—differentiable function 
at t, where t > 0. Then the conformable fractional derivative of f of order a 
is defined as 


fa]— fal—a)) _ ¢(fa]— 
‘pg iat oe ee (3) 


e0 E 


where [a] is the smallest integer greater than or equal to a. 
When the conformable fractional derivative of f of order a exists, we say 
f is a—differentiable and we write f‘© (t) for Ty(f)(t). 


Remark 1. [18] As a consequence of Definition 1, one can easily show 
that 


Tol f)(t) = Foto) pls (), (4) 
where a € (n,n + 1], and f is (n + 1)—differentiable at t > 0. 


Theorem 1. [18] Let a € (0,1], and f,g be a—differentiable at a point 
t > 0. Then 
1. To(af + bg) =aTa(f)+bTa(g), for all a,be R, 
2. Ta(fg) = fTa(g) + 9Ta(f). 

In [1], Abdeljawad was defined the left and right conformable fractional 
derivative. Since the left fractional derivative on [0,0o) is the conformable 
fractional derivative, we can have the following theorems according to [1]. 


Theorem 2. (Chain Rule) Assume f,g : (0,00) —> R be a—differentiable 
functions, where 0 < a < 1. Let A(t) = f(g(t)). Then A(t) is a—differentiable 
and for all t with t 4 0 and g(t) 4 0 we have 


40 E. Hesameddini and E. Asadollahifard 


(Tah)(t) = (Taf)(9(t))(Tag) (t)g(t)?™. 
If t = 0, then 


(T.h)(0) = jim (Taf)(9(t) Fagg)? 


Theorem 3. Let f : (0,00) — R be twice differentiable on (0,00) and 
0<a,8<1such that 1<a+ $< 2. Then 


(ToT f)(t) = Tore f(t) + (1 — 8)tTaf(t). 


2.2. Sinc function 
The sine function is defined on the whole real line, —co < x < cw, by 


sinc(x) = { a 2 “a 7 


For each integer & and the mesh size h, the translated sinc basis function is 
defined as 


x—kh 


8(k, h)(a) = sinc( ). 


If a function f(a) is defined on the real axis, then for any h > 0, the Whittaker 
cardinal expansion of f(a) is as follows 


a—kh 


c(f, h)(x) = ys f (kh) sinc( 


k=—0o 


), 


whenever this series converges. The properties of Whittaker cardinal expan- 
sion are derived in the infinite strip D, of the complex w-planes where for 
d>0 


Dy = {w =t+is:|s|<d< 5}. 


These properties have been studied thoroughly in [24]. In order to approxi- 
mate on the finite interval (a,b), which is used in this paper, we consider the 


one-to-one conformal map w = ¢(z) = In(4=$), which maps the eye-shaped 
domain 
z-a T 
De={z= Ly : d<s= 
p={z=axut+iy largp— | < Soi 


onto the infinite strip D,. The basis functions on (a,b) are taken to be the 
composite translated sinc functions 
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o(a) — kh 
h 


$4(x) = 8(k, h)od(x) = sinc( ), kEZ (5) 


where s(k, h)od(x)is defined by s(k, h)((x)). 
Let / = ¢-|. We define the range of ~ on the real line as 


T= {w(w) € De: —co < w < co}. 


For the uniform grid {kh}? on the real line, the image which corresponds 
to these nodes is denoted by 


1 


tp = (kh) = ee k= Dict (6) 


For discretizing the problem we need the following definition and theorems. 


Definition 2. [24] Let Lg(Dz) be the set of all analytic functions, for 
which there exist a constant, C, such that 


ela 
BS OT Tote 


2eDn, UA pS 


where p(z) = e%) . 
Theorem 4. [21] Let y € Lg(Dz), N be a positive integer and h be selected 
by the formula 


rd 


b= (ae) (7) 


then there exists a positive constant c,, independent of N, such that 


N 


1 
supzerly(z)— > y(z)sG, h)og(2)| < cre ON)”. 
j=—-N 


Theorem 5. [21] Let ¢ be a conformal one-to-one map of the simply con- 
nected domain Dg onto Dg.Then 


i_ 
PanSles= (4402 


0 k=J, 
—1)G-*) , 
eye hay 


a | 


kj = Fp loe(w)llens, = 
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Pp 1 Lee k ' 
=); 

5(2) = 7 |Sk(e = = WZ 2 eget 

kj [ (x)| j eae kj. 


3 Method of Solution 


Consider equation (1) in J = [0,1] where D%y denotes the fractional deriva- 
tive which is defined in (3) ie. DOMy = y, 

The approximate solution of equation (1) based on the sinc basis functions 
(5), should satisfy the initial conditions (2). But this basis functions do not 
have a derivative when x tends to 0 or 1 so we modify them as 


w(x)sx(x), (8) 


where w(x) = (a(1 — ))&"-» [6]. 
In order to approximate the solution, we construct a polynomial p(x) that 
satisfies initial conditions [6]. So the approximate solution is represented by 


yn (x) = un(x) + p(2), (9) 
where 
N 
un(x)= S~ cyw(x)sx(2), (10) 
k=—N 
and 
p(x) = ap taya+...+Am2z™, m-1L<v<m. (11) 


The unknown coefficients ao, @1,..., @m and {c,}/__ y are determined by sub- 
stituting yy (x) into equation (1) and evaluating the result at the sinc points 


eh 
~ 1l+eih 


x; j=—-N-1,...,N. (12) 


Notice that according to Theorem 1 and Remark 1, we have 


(w(x) sx(ax)) = attlel-% (q(x) sy (2))OF), n<a<ntl, (13) 
so 
uy (x) = BRL vee (w(x) se(x)). (14) 


Also it should be noted that when x tends to 1 or 0, we have 
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uns) = aya) =.= ul?) (q) = 0. 


Using equations (13) and (14), one can obtain 


y) (x,;) — u© (a;) + pO (x;), j=—-N-1,...,N. (15) 


Now by substituting this equation into equation (1), we obtain the following 
system of algebraic equations which can be solved for unknowns 
un (23) = F(e3, yn (23), YN” (@), 9" (Bs), —N-1LSGSN, 


(0) =di, §=0,1,.m—1. 


4 Applications and results 


In this section, we solve some examples by the presented method and com- 
pare the numerical results with the exact solutions and some earlier works. 


Example 1. As the first example, we consider the following nonlinear frac- 
tional initial value problem [5] on [0, 1] 


y" (x) +y??(a) +y?(x) =24, (0) =y'(0) =0, y"(0)=2, (16) 


whose exact solution is y(2) = x?. Following the procedure of the presented 
method, we consider the following approximate solution 


k=—N 


where w(x) = 2?(1— 2x). By substituting this approximate solution into 
equation (16) and evaluating at sinc points (12), we arrive at the following 
nonlinear system of algebraic equations which can be solved for unknown 
coefficients 


Gaz + OR yen{wld + 6 (Bwid) + 3wib! + wh") + 5. (3wi(4)?4+ 
30 b4b4) + bgp Wy (O))>} + 6a} ag + DE_wenfap Pw! Jy) + wy (OY 5., + 
3616552 + (#i)98)} + (a0 + a12; + 0203 + aga? + SP_yexw;d)? = 0 
j=—-N-1,...,N, 
y(0)=0Sa,=0, y’'(0)=0Sa,=0, y"(0)=2> a2 =1. 
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According to relation (7), by taking d = 5 and 6 = 2, we have h = Wink 
Then by applying the well known Newton method with starting points cz, = 
0,k = -N,...,N, agp = a1 = az = 0,a2 = 1, we obtain c, = 0,k = —N,...,N 
and aj = a1 = a3 = 0,a2 = 1. So the approximate solution is yy(x) = 2?, 


which is the exact solution. 


Example 2. Consider the fractional Ricatti equation on [0,1] 
yO (x) = 2y(a) - y?(@) +1, O<aK<1, y(0) =0. 


For a = 1, the exact solution of this equation is y(t) = 1+ V2tanh(V/2t + 
sln( vi-t)). Consider the following approximate solution based on the sinc 


collocation method 


yn (x) = SR__ yeesy(x) + ao + aye. 


Odibat and Momani [20], solve this equation by using the modified homotopy 
perturbation method. Also in [4], this equation is solved by the Chebyshev 
wavelet operational matrices of fractional integration. For comparison, the 
results of this method are presented in Tables 1 and 2 with 192-set of Block 
Pulse Functions ( Chebyshev wavelets was expanded into an 192-term block 
pulse functions). 

In Table 1, the results of the presented method with N = 1 for a = 0.5 and 


Table 1: Numerical results with comparison to [4, 20] for a = 0.5 and a = 0.75 in 
Example 2 
a=0.5 a = 0.75 

x | Ours[N=1] [4] [20] Ours[N=1] [4] [20] 
0.1 | 0.3956920 0.592756 | 0.321730 | 0.2321153 | 0.310732 | 0.216866 
0.2 | 0.9184524 | 0.9331796 | 0.629666 | 0.4961556 | 0.584807 | 0.428892 
0.3 | 1.2973611 | 1.1739836 | 0.940941 | 0.7523005 | 0.822173 | 0.654614 
0.4 | 1.5802323 | 1.38466546 | 1.250737 | 0.9998683 | 1.024974 | 0.891404 
0.5 | 1.7987123 | 1.4738876 | 1.549439 | 1.2372036 | 1.198621 | 1.132764 
0.6 | 1.9690794 | 1.5705716 | 1.825456 | 1.4604023 | 1.349150 | 1.870240 
0.7 | 2.0982657 1.646199 | 2.066523 | 1.6619744 | 1.481449 | 1.594278 
0.8 | 2.1867519 1.706880 | 2.260633 | 1.8278045 | 1.599235 | 1.794879 
0.9 | 2.2352250 1.756644 | 2.396839 | 1.9347648 | 1.705303 | 1.962239 
1.0 | 2.8926026 1.798220 | 2.466004 | 2.0825668 | 1.801763 | 2.087384 


a = 0.75 are compared with earlier works [4, 20]. We see that our results are 
in a good agreement with them. For a = 1, the results are presented in Table 
2. It is clear that by increasing N, the approximate solution becomes more 
and more accurate and for N = 85 the exact solution is obtained whereas 
Refs [4, 20] can not reach the exact solution. In Figure 1. the approximate 
solution for different values of a is shown. Numerical results show that as 
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Table 2: Numerical results with comparison to [4, 20] for a = 1 in Example 2 

x | Ours[N=10] | Ours[N=50] | Ours[N=80] [4] [20] Exact 
0.1 0.1134865 0.1103047 0.110295 0.1103111 | 0.110294 | 0.110295 
0.2 | 0.2458331 0.2419881 0.241977 0.241995 | 0.241965 | 0.241977 
0.3 | 0.3993884 0.3951178 0.395105 0.395123 | 0.395106 | 0.395105 
0.4 | 0.5726929 0.5678265 0.567812 0.567829 | 0.568115 | 0.567812 
0.5 | 0.7610790 0.7560297 0.756014 0.756029 | 0.757564 | 0.756014 
0.6 | 0.9589295 0.9535820 0.953566 0.953576 | 0.958259 | 0.953566 
0.7 1.1581332 1.1529646 1.152949 1.152955 | 1.163459 | 1.152949 
0.8 | 1.3514117 1.3463785 1.366364 1.346365 | 1.365240 | 1.366364 
0.9 1.5314497 1.5269249 1.526911 1.526909 | 1.554960 | 1.526911 

1 1.6949935 1.6895135 1.689498 1.689494 | 1.723810 | 1.689498 


@ approaches to its integer value, the solution of fractional order differential 


equation approaches to the solution of integer order differential equation. 


Figure 1: Approximate solution of Example 2 for different values of a 


Example 3. [3] As the last example, consider the following inhomogeneous 
Bagley-Torvik equation 


subject to initial conditions 


y" (a) +y%) (2) + y(z) =1+2, 


The exact solution of this equation is y(v) = 1+ 2. 
In a same manner of last examples, by considering the approximate solution 


as 


yn (x) = SR__ yegw(ax)sp(x) + ap + a2 + ax, 
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where w(x) = x(1 — 2x), one can obtain yy(x) = 1+ which is the exact 
solution. 


5 Conclusion 


In this work the sinc-collocation method is used to approximate the solution 
of multi-order fractional differential equations with initial conditions. This 
method converts the FDEs into a system of algebraic equations which can be 
solved more easier. In this work, the fractional derivatives are described in 
the sense of new definition which makes us able to solve fractional differential 
equation directly by the sinc method for the first time. Also this method can 
be applied to other types of FDEs easily. Several examples are included to 
demonstrate the reliability and efficiency of our method. 
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